function B = EmCir3D(x_d,y_d,z_d,l,ai,ao)

r_d=sqrt(x_d.^2+y_d.^2);
if r_d>0
    cos_theta_d=x_d/r_d;
    sin_theta_d=y_d/r_d;
    
    B_R=triplequad(@(r,z,theta)EMCirFuncR...
        (r,z,theta,r_d,z_d),ai/2,ao/2,-l/2,l/2,0,2*pi);
else
    B_R=0;
    cos_theta_d=0;
    sin_theta_d=0;
end

B_Z=triplequad(@(r,z,theta)EMCirFuncZ...
    (r,z,theta,r_d,z_d),ai/2,ao/2,-l/2,l/2,0,2*pi);

B=[B_R*cos_theta_d B_R*sin_theta_d B_Z];
end
function B = EMCirFuncR(r,z,theta,r_d,z_d)

nu0 =  4*pi*10^-7;
R=sqrt(r.^2+r_d.^2-2*r.*r_d.*cos(theta)+(z-z_d).^2);
Nom=(z_d-z).*cos(theta).*r;
Denom=R.^3;
Const=nu0/(4*pi);

B=Const*Nom./Denom;
end

function B=EMCirFuncZ(r,z,theta,r_d,z_d)

nu0 =  4*pi*10^-7;
R=sqrt(r.^2+r_d.^2-2*r.*r_d.*cos(theta)+(z-z_d).^2);
Nom=r.*(r-r_d.*cos(theta));
Denom=R.^3;
Const=nu0/(4*pi);
B=Const*Nom./Denom;
end
